{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e16fe9fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1aa816c2a30>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAu90lEQVR4nO3deXhU5fn/8fednYSEJSQhJIQAWSCAbBHBBWUVUItaa8Vq8ddWagXFpbZYrbXW1q3aqkUUAaVuuCsqCgqoLLIEZIeQEAgJWwKBkBCyTPL8/pjBbxoTkjCTnFnu13XNlTnbnPtAmA/PWZ5HjDEopZTyXX5WF6CUUspaGgRKKeXjNAiUUsrHaRAopZSP0yBQSikfF2B1AeeiU6dOJjEx0eoylFLKo2zYsOGoMSaq7nyPDILExEQyMjKsLkMppTyKiOTWN19PDSmllI/TIFBKKR+nQaCUUj5Og0AppXycBoFSSvk4lwSBiMwTkQIR2dbAchGR50QkW0S2iMigWsvGiUimY9kMV9SjlFKq6VzVIngVGHeW5eOBZMdrCjALQET8gZmO5WnAJBFJc1FNSimlmsAlzxEYY74VkcSzrDIR+K+x93m9RkTai0gskAhkG2NyAERkgWPdHa6oSzVPeVU1OYWnyCoooehUJRW2GiqqaggMEKLDQ4gODyYxMoyuHdsgIlaXq5RykdZ6oCwOyKs1ne+YV9/8C+r7ABGZgr01QUJCQstU6WMqbTVk5BbxdWYhX2cWkFVQSlOGp+jUNoiBCR0Y2iOSK8+LJSYipOWLVUq1mNYKgvr++2jOMv/HM42ZDcwGSE9P19F0nJBXVMbra3J5OyOPE2VVBPoLF3SPZHzfWJJj2pIU3Zbo8BBCAv0I8vejwlZDYUkFBSUVZBeUsiH3OBtyi/hyxxEe/WwHF/aM5NqB8VzVvwtBAXr/gVKeprWCIB/oWms6HjgIBDUwX7WAnYdO8vSS3SzddQQ/EcamxXD1wDguSupE2+CGfxUC/P0ICw4gsVMYQ7p35MYL7C2yPYWlfLzpIB9vOsC9727m6SWZTBnegxuGJBAS6N9ah6WUcpK4aqhKxzWCT40xfetZdgUwDZiA/dTPc8aYISISAOwGRgEHgPXAjcaY7WfbV3p6utG+hpruwInTPL0kkw+/P0BESCCTh3Vj0gUJxLZr45LPN8bwbdZRZi7LZt2+IqLCg7l/fC+uGRin1xKUciMissEYk153vktaBCLyFnAZ0ElE8oG/AIEAxpgXgUXYQyAbKAP+n2OZTUSmAYsBf2BeYyGgmq66xvDKqr08tTgTA0y5pAe3X5ZEu9BAl+5HRLg0JYpLU6JYm3OMxz7fxT3vbGbB+jz+NrEvqZ3DXbo/pZRruaxF0Jq0RdC4nMJS7ntvCxtyjzO6dzR/ndiXuPauaQE0pqbG8HZGHk98sYvSchu/vzyVKZf0wM9PWwdKWalFWwTKvby/IZ8HPtpKcIA///p5f64e0LqnaPz8hElDEri8T2ce/Ggrj3++i1XZR3n6+v5Eh+sdRkq5G73Fw4tU2mp46ONt3PvuZgZ0bc+Su4dzzcB4y87TdwwLYuaNg3js2n6s31fEhGdXsG5vkSW1KKUapkHgJY6WVjDp5TX897tcbr2kO6//+gK3uL9fxN46+GTaxUSEBHLTnLV8sDHf6rKUUrVoEHiBvKIyrpu1mu0Hi3l+0kAeuCKNAH/3+qtNjgnnw9svYnC3Dtzzzmb+uTiTmhrPuz6llDdyr28L1WzbDxZz7azVHC+r4o3fDOWq/l2sLqlB7UIDmf+rIfw8vSv/WZ7Nnz7cSrWGgVKW04vFHmxDbhG3zFtP25AA3rxtGMkx7n+bZlCAH4//tB/REcE8vyyb8qpq/vmz/m7XglHKl2gQeKiN+48zed56osKDeeM3F9CllW4NdQUR4d6xqYQE+vPU4kwqbDU8e8NA7Z5CKYvovzwPtCX/BJPnriOybRBv3TrUo0KgtqkjkvjzlWl8vu0wd7+zSU8TKWURbRF4mB0HT3LTnLW0Cw3kzVuH0rmd9XcGOePXF3enpsbw90U7iQgJ5B/X9NVuKZRqZRoEHiT/eBm3vLKOsOAA3rp1aKs9KdzSbh3egxOnK5m5fA/tQwP547heVpeklE/RIPAQxWVV3PLKek5XVfPebRfStWOo1SW51O/HpnKirIpZX++hU9tgfn1xd6tLUspnaBB4gPKqam79bwb7j5Ux/1dDvLITNxHhkYl9KTpVyaOf7SAxMpRRvWOsLkspn6AXi92cMYb7P9jKun1F/PP6/gzrGWl1SS3G30945voB9O3Sjjve+p4dB09aXZJSPkGDwM3NXbmXD78/wL1jUviJGz8s5iptgvyZMzmdiJBAfjN/PQUny60uSSmvp0HgxlZkFfKPRTsZ37cz00YmWV1Oq4mJCGHO5HSOl1Vx+xsbqaqusbokpbyaBoGbyj12imlvfk9KTDj//Fl/n7ulsm9cO5647jwyco/z2KJdVpejlFdzSRCIyDgRyRSRbBGZUc/y+0Rkk+O1TUSqRaSjY9k+EdnqWKajzWC/OPy71zcCMPvmdMLOMp6wN/tJ/y7ccmEi81bt5dMtOpS1Ui3F6SAQEX9gJjAeSAMmiUha7XWMMU8ZYwYYYwYA9wPfGGNqd0w/wrH8RyPn+KJHP9vBjkMneeb6/iREetdtos31pwm9GdytA394bwtZR0qsLkcpr+SKFsEQINsYk2OMqQQWABPPsv4k4C0X7NcrfbrlIK+v2c+U4T309knsndS98ItBhAb5M/XNjZRXVVtdklJexxVBEAfk1ZrOd8z7EREJBcYB79eabYAlIrJBRKY0tBMRmSIiGSKSUVhY6IKy3c++o6eY8f5WBia0577LU60ux23ERITwzPUD2H2klL9/ttPqcpTyOq4IgvquYjbUe9hVwKo6p4UuMsYMwn5qaaqIDK9vQ2PMbGNMujEmPSoqyrmK3VBVdQ3TF3yPv5/wnxsHEajdMv+P4SlR3HpJd15bk8uS7YetLkcpr+KKb5t8oGut6XigoSt7N1DntJAx5qDjZwHwIfZTTT7n+WXZbM4v5rFr+3lNH0Ku9vvLU+nTJYI/vL+FQ8WnrS5HKa/hiiBYDySLSHcRCcL+Zb+w7koi0g64FPi41rwwEQk/8x4YC2xzQU0eZeP+48xcns21g+KY0C/W6nLcVnCAP89NGkhFVQ2/f3ezDnWplIs4HQTGGBswDVgM7ATeMcZsF5HbROS2WqteAywxxpyqNS8GWCkim4F1wGfGmC+crcmTnKqwcffbm+gcEcLDP+ljdTlur2dUWx68sjerso/xxtpcq8tRyiu45AZ1Y8wiYFGdeS/WmX4VeLXOvBygvytq8FR/X7ST/UVlvD1lGBEhgVaX4xFuHJLAF9sO849Fu7gkOYrETmFWl6SUR9MrkhZalX2UN9fu59ZLejCke0ery/EYIsKT151HgL9w33ubdWQzpZykQWCRUxU2/vj+Fnp0CuOeMSlWl+NxYtu14eGr+rB+33FeWbXX6nKU8mgaBBZ5anEmB06c5snrziMk0N/qcjzStYPiGN07hn8uyST32KnGN1BK1UuDwALr9hbx6up9TB6WSHqinhI6VyLCo1f3JdDPjz99uBVj9BSRUudCg6CVVdiqmfH+Frp2bMMfxunTw87q3C6EP47vxarsY7y3Id/qcpTySBoErWzW13vIOXqKv1/dj9Ag3+xV1NVuHJLA+YkdePSznRSWVFhdjlIeR4OgFe0pLOWF5Xv4Sf8uDE/xvm4yrOLnJzx27Xmcrqzmr59st7ocpTyOBkErMcbwwIdbCQn048Ere1tdjtdJim7L1BFJfLrlECuyvLNTQqVaigZBK3l/4wHW5BQxY3xvosNDrC7HK/320h4kRoby0MfbqbBpd9VKNZUGQSsoLqviH4t2MiihPTec37XxDdQ5CQn055GJfdl79BSzv8mxuhylPIYGQSt45stMTpRV8ujV/fDz862xh1vb8JQorugXy3+WZ5NXVGZ1OUp5BA2CFrbj4EleW5PLzUO7kdYlwupyfMKfr0wjwE94eKFeOFaqKTQIWpAxhr8s3Eb70CDuGaPPDLSWzu1CmD46maW7Cli+q8DqcpRyexoELeijTQdYv+84fxyXSrtQ7Vm0Nd1yYXd6dArjb5/uoNJWY3U5Srk1DYIWUlph4x+LdtG/a3t+NlgvELe2oAA//nxVGjlHT/Hqau2UTqmzcUkQiMg4EckUkWwRmVHP8stEpFhENjleDzV1W0816+tsCksqePiqNL1AbJERqdGM7BXNc0uzKSgpt7ocpdyW00EgIv7ATOyDz6cBk0QkrZ5VVxhjBjhejzRzW4+SV1TGyyv2cs3AOAYmdLC6HJ/25yvTqLBV89QXmVaXopTbckWLYAiQbYzJMcZUAguAia2wrdt6/PNd+Itop3JuoHunMH51UXfe3ZDPtgPFVpejlFtyRRDEAXm1pvMd8+oaJiKbReRzETkzOG9Tt/UY6/YW8dnWQ9x2aU9i27WxuhwFTB2ZRMewIB79bId2Va1UPVwRBPWdAK/7r20j0M0Y0x94HvioGdvaVxSZIiIZIpJRWOiefcnU1Bge+XQ7se1CmDK8h9XlKIeIkEDuHp3MmpwivtxxxOpylHI7rgiCfKD2bTHxwMHaKxhjThpjSh3vFwGBItKpKdvW+ozZxph0Y0x6VJR79tz50aYDbDtwkj+MS6VNkI465k4mDUkgKbotj32+S28nVaoOVwTBeiBZRLqLSBBwA7Cw9goi0llExPF+iGO/x5qyracor6rmn4sz6RsXwcT+Hn12yysF+PvxwITe7D16ijfW5lpdjlJuxekgMMbYgGnAYmAn8I4xZruI3CYitzlWuw7YJiKbgeeAG4xdvds6W5MVXlm1j4PF5fxpQm+9XdRNXZYaxcVJnfj3V1kUl1VZXY5SbkM88eJZenq6ycjIsLqMHxSdquTSJ5czpHtH5t5yvtXlqLPYcfAkVzy/ginDe3D/eB0XQvkWEdlgjEmvO1+fLHaB55ZmcarSxozxvawuRTUirUsE1wyIs7fgTpy2uhyl3IIGgZNyj9nPOf/8/ASSY8KtLkc1wT1jU8DAM1/utroUpdyCBoGTnl6ymwA/P+4enWx1KaqJ4juEMvnCbry/MZ9dh09aXY5SltMgcMK2A8Us3HyQX12cSHSEDj/pSaaOSCI8OIAntesJpTQInPHk4kzahwby20t7Wl2Kaqb2oUHcPiKJZbsKWJtzzOpylLKUBsE5Wr3nKN/uLmTqZUlEhOhYA57olgsTiYkI5snFmdr1hPJpGgTnwBjDE19k0qVdCDcP62Z1OeochQT6M31UChtyj7NMRzJTPkyD4Bws3n6EzXknuGt0CiGB2pWEJ/tZejyJkaE8tTiTmhptFSjfpEHQTNU1hme+zKRHVBjXDtKuJDxdoL8f94xNZdfhEj7ZUm83V0p5PQ2CZlq4+QC7j5Ry75hUAvz1j88bXNkvlrTYCJ5esls7pFM+Sb/JmqGquoZ/fZlFWmwE4/t2troc5SJ+fsJ9l6eyv6iM9zbkW12OUq1Og6AZ3snIY39RGfddnqody3mZy1KjGNytA88vy6K8qtrqcpRqVRoETVReVc1zS7MY3K0Dl6W653gI6tyJCPeOTeFQcTlvrt1vdTlKtSoNgiZ6fU0uR05W8PuxqTiGVlBe5sKenbiwZyQvfJ1NWaXN6nKUajUaBE1QVmnjxW/2cFFSJMN6RlpdjmpB945N4WhpJfNX6+A1yndoEDTBf7/L5WhpJfeMSbG6FNXCBnfryIjUKF78Zg8ny3XwGuUbXBIEIjJORDJFJFtEZtSz/BcissXxWi0i/Wst2yciW0Vkk4i4z2gzDqUVNl76Zg+XpkQxuFtHq8tRreCeMakUn67ilZX7rC5FqVbhdBCIiD8wExgPpAGTRCStzmp7gUuNMecBfwNm11k+whgzoL6Rc6z2ysq9HC+r0taAD+kX344xaTHMWZlD8WltFSjv54oWwRAg2xiTY4ypBBYAE2uvYIxZbYw57phcA8S7YL8trvh0FS+vyGF07xj6d21vdTmqFd01OpmSchtzV+61uhSlWpwrgiAOyKs1ne+Y15BfA5/XmjbAEhHZICJTGtpIRKaISIaIZBQWFjpVcFPNW7mXk+U27h6jg874mj5d2jGuT2deWbmXE2WVVpejVItyRRDUdy9lvb13icgI7EHwx1qzLzLGDMJ+ammqiAyvb1tjzGxjTLoxJj0qquXv4y8uq2Leyr2M69OZPl3atfj+lPu5a0wyJRU25qzQVoHybq4Ignyga63peOBHvXeJyHnAHGCiMeaHkUCMMQcdPwuAD7GfarLc3JU5lFTYmK5DUPqsXp0juKJfLK+s2kvRKW0VKO/liiBYDySLSHcRCQJuABbWXkFEEoAPgJuNMbtrzQ8TkfAz74GxwDYX1OSUE2WVvLJqH+P7dqZ3bITV5SgLTR+dTFlVNXNW5FhdilItxukgMMbYgGnAYmAn8I4xZruI3CYitzlWewiIBF6oc5toDLBSRDYD64DPjDFfOFuTs+au3EtJhY07R2lrwNelxIRzRb9Y5q/ep60C5bUCXPEhxphFwKI6816s9f43wG/q2S4H6F93vpXOtAYm9NPWgLK7c1Qyn209xJwVOfxhXC+ry1HK5fTJ4jrmrtxLqbYGVC21WwXHtVWgvJAGQS21WwO9OmtrQP2fO0c5rhWs1GsFyvtoENQyT1sDqgFnWgWvrtJWgfI+GgQOxWVVvLJqH+P6aGtA1U9bBcpbaRA4zFuldwqps0uJCWdC31jmr87Vp42VV9EgwN6n0LxVexmbFkNaF20NqIbdMSqJ0gob87QPIuVFNAiA+av3UVKurQHVuF6dI+x9EK3apz2TKq/h80FQUl7F3JV7Gd07hr5x2qeQatwdo5IoqbDxyiptFSjv4PNB8N/vcik+XcWdo5KsLkV5iD5d7OMV2Hun1VaB8nw+HQSnKmzMWZHDiNQozotvb3U5yoNMH5XMyXIb81fts7oUpZzm00Hw+ppcjpdV6bUB1Wx949oxqlc0c1fZnz1RypP5bBCcrqxm9rc5DE+JYmBCB6vLUR7ojlHJnCir4rXvcq0uRSmn+GwQvLE2l2OnKpmu1wbUORrQtT2XpkQxZ0UOZZXaKlCeyyeDoLyqmpe+zeHCnpEM7tbR6nKUB7tzVDLHTlXy5tr9Vpei1DnzySB4e30ehSUVem1AOW1wtw5clBTJi9/kUF5VbXU5Sp0TnwuCCls1s77ew5DEjgztEWl1OcoL3DkymaOlFby1TlsFyjO5JAhEZJyIZIpItojMqGe5iMhzjuVbRGRQU7d1tXcz8jl8slxbA8plLugRyQXdO/LSNzlU2LRVoDyP00EgIv7ATGA8kAZMEpG0OquNB5IdrynArGZs6zKVthpmfb2HgQntuShJWwPKde4clczhk+W8m5FvdSlKNZsrWgRDgGxjTI4xphJYAEyss85E4L/Gbg3QXkRim7ity3z4fT4HTpzmzlHJiEhL7Ub5IPuNBx2Y9fUeKm01VpejvNDxU5XcPHct2w4Uu/yzXREEcUBerel8x7ymrNOUbQEQkSkikiEiGYWFhedUaGFJBendOnBZStQ5ba9UQ0SEO0YmceDEaT78XlsFyvXmrdrLiqyjBAW4/tKuKz6xvv9amyau05Rt7TONmW2MSTfGpEdFndsX+bSRybz922HaGlAt4tKUKM6Lb8fM5XuwVWurQLlO8ekqXnUMo5sSE+7yz3dFEOQDXWtNxwMHm7hOU7Z1KX8/DQHVMkSEO0cms7+ojI82teivsfIxr67aR0mFjWkjWuYmF1cEwXogWUS6i0gQcAOwsM46C4FfOu4eGgoUG2MONXFbpTzGqN7RpMVGMHN5NtU19TZulWoWe1f5OYxpwYGznA4CY4wNmAYsBnYC7xhjtovIbSJym2O1RUAOkA28DNx+tm2drUkpq4gId45KYu/RU3y6RVsFynn//S6Xk+U27hzZcre8B7jiQ4wxi7B/2dee92Kt9waY2tRtlfJkY9M6kxoTzvPLsrnqvC746elIdY5qd5XfL77lBs7yuSeLlWppfn7CHaOSyC4o5fNth60uR3mwN9bau8q/o4UfgNUgUKoFjO8bS8+oMJ5flkWNXitQ5+BMV/mXJHdiUAt3la9BoFQL8PcT7hiZzK7DJSzZccTqcpQHenPdfo6WVjK9FbrD0SBQqoVceV4s3TuF8dzSLOyXyZRqmvKqal78Zg8X9owkPbHlu8rXIFCqhQT4+zF1RBI7Dp1k6c4Cq8tRHqS1u8rXIFCqBV09oAsJHUN5VlsFqol+6Cq/e+t1la9BoFQLCvD3Y9qIJLYeKObrzHPrI0v5lnfOdJXfgs8N1KVBoFQLu2ZQHPEd2vBvbRWoRlTYqpm1PPuHke9aiwaBUi0s0HGtYHPeCb7Zra0C1bD3NuRzsLic6a3cVb4GgVKt4KeD4olr30avFagGVdpqeGG5feCsS5I7teq+NQiUagVBAX787rKefL//BCuyjlpdjnJD722wD5zV2q0B0CBQqtX8LD2eLu1C+PdXu7VVoP5Hpa2Gmcuz6d+1PZdaMHCWBoFSrSQ4wJ/bRySxUVsFqo4PNtpbA3dZNIyuBoFSrehMq0CvFagzKm01/Gd5Nv3j23FZqjXD6GoQKNWKzrQKNuQeZ2W2tgoUvL8xn/zjp7lrTIplw+hqECjVys60Cv71pV4r8HWVthr+syybAV3bc5kF1wbOcCoIRKSjiHwpIlmOnz/qK1VEuorIchHZKSLbRWR6rWUPi8gBEdnkeE1wph6lPEHtawXf6rUCn3bmTqG7RltzbeAMZ1sEM4ClxphkYKljui4bcK8xpjcwFJgqImm1lv/LGDPA8dKRypRPuD69K3Ht22irwIeduVNoYII1dwrV5mwQTATmO97PB66uu4Ix5pAxZqPjfQn2sYnjnNyvUh4tKMCPaSOT2JR3Qvsg8lHvbshztAasuzZwhrNBEGOMOQT2L3wg+mwri0giMBBYW2v2NBHZIiLz6ju1VGvbKSKSISIZhYX6D0d5vusGxxPfoQ3/0ucKfE55VTX/WZbNoIT2DG/lp4jr02gQiMhXIrKtntfE5uxIRNoC7wN3GWNOOmbPAnoCA4BDwNMNbW+MmW2MSTfGpEdFWduMUsoVAv39uHNkMlvyi3W8Ah+zYN1+DhWXc+/YVMtbA9CEIDDGjDbG9K3n9TFwRERiARw/6/1tFpFA7CHwhjHmg1qffcQYU22MqQFeBoa44qCU8hTXDIqjW2Sotgp8yOnKamZ+vYcLunfkwp6t18Po2Th7amghMNnxfjLwcd0VxB53c4Gdxphn6iyLrTV5DbDNyXqU8ihnWgXbD55k8fbDVpejWsHra3IpLKngHgufG6jL2SB4HBgjIlnAGMc0ItJFRM7cAXQRcDMwsp7bRJ8Uka0isgUYAdztZD1KeZyrB8bRMyqMZ77cTXWNtgq82akKGy9+s4eLkzpxQSuNPtYUAc5sbIw5BoyqZ/5BYILj/Uqg3tgzxtzszP6V8gb+fsLdY1KY9ub3fLL5IFcP1JvqvNX87/Zx7FQld49JsbqU/6FPFivlBib0jaVX53D+/dVuqqprrC5HtYDi01W89E0OI1KjGNytwRskLaFBoJQb8PMT7h2byr5jZXywMd/qclQLmLMih+LTVdw7NtXqUn5Eg0ApNzG6dzT9u7bnuaXZVNiqrS5HudDR0grmrtzLFf1i6RvXzupyfkSDQCk3ISL8fmwKB06c5s21+60uR7nQrK/3UF5V7XbXBs7QIFDKjVyc1IlhPSL5z7JsSitsVpejXOBQ8WleW5PLtYPiSYpua3U59dIgUMqNiAj3jUvl2KlK5q3ca3U5ygWeW5qNMYbpo5KtLqVBGgRKuZlBCR0YkxbDy9/mcPxUpdXlKCfsKSzlnYw8bhySQNeOoVaX0yANAqXc0O/HplJaaWPWN3usLkU54eklmQQH+DFtpPu2BkCDQCm3lNo5nGsGxjF/9T4OnjhtdTnqHGzOO8GirYf5zSU9iAoPtrqcs9IgUMpN3T06BWPgX1/utroU1UzGGJ74Yhcdw4K49ZLuVpfTKA0CpdxU146h3DysG+9vzCfzcInV5ahmWJF1lNV7jjFtRBLhIYFWl9MoDQKl3Ni0EUmEBQfwxBe7rC5FNVFNjb01ENe+Db8YmmB1OU2iQaCUG+sQFsTvLuvJsl0FrMk5ZnU5qgk+2nSA7QdP8odxqQQH+FtdTpNoECjl5n51UXc6R4Tw2Oe7dPAaN1deVc0/F2fSL64dV53XxepymkyDQCk3FxLozz1jUticd4JPtxyyuhx1FvNW7eVgcTl/mtAbPz/3GHSmKZwKAhHpKCJfikiW42e9fauKyD7HADSbRCSjudsr5et+OjieXp3DeeKLXZRXaYd07uhYaQWzlu9hdO9ohrnJEJRN5WyLYAaw1BiTDCx1TDdkhDFmgDEm/Ry3V8pn+fsJD16RRv7x07y6ep/V5ah6PL8sm7KqamaM72V1Kc3mbBBMBOY73s8Hrm7l7ZXyGRcnd2Jkr2hmLsvmWGmF1eWoWrILSnhtTS4/P78rSdHhVpfTbM4GQYwx5hCA42d0A+sZYImIbBCRKeewvVIK+NOEXpRVVfPvr7KsLkXV8vfPdhLquJbjiRods1hEvgI617PogWbs5yJjzEERiQa+FJFdxphvm7E9jgCZApCQ4Bn35irlaknR4dw4JIE31+3nl8O6kRzjef/79DZfZxawPLOQByb0plNb9+5KoiGNtgiMMaONMX3reX0MHBGRWADHz4IGPuOg42cB8CEwxLGoSds7tp1tjEk3xqRHRUU15xiV8ip3jU4mLMifRz7dobeTWsxWXcOjn+0kMTKUyRcmWl3OOXP21NBCYLLj/WTg47oriEiYiISfeQ+MBbY1dXul1P+KbBvM3WNSWJF1lK92Nvh/J9UK3ly3n+yCUv40oTdBAZ57N76zlT8OjBGRLGCMYxoR6SIiixzrxAArRWQzsA74zBjzxdm2V0qd3U1Du5Ec3Za/fbpDbye1yPFTlTzz5W4u7BnJmLQYq8txSqPXCM7GGHMMGFXP/IPABMf7HKB/c7ZXSp1doL8ff7mqDzfNXcvclXuZOiLJ6pJ8zlNLMikpt/GXq/og4jkPj9XHc9sySvm4i5M7cXmfGGYuz+ZwcbnV5fiULfkneGvdfiYPSyS1s+dfsNcgUMqDPXhFGtU1hr8v2ml1KT6jpsbw0MfbiQwL5q4x7j3yWFNpECjlwbp2DGXqiCQ+2XyQlVlHrS7HJ7y3MZ9NeSe4f3wvIjxgrIGm0CBQysNNGd6DxMhQHvp4GxU2vXDckk6UVfLE57sY3K0D1wyMs7ocl9EgUMrDhQT688jEvuQcPcXL3+ZYXY5Xe/zzXZw4XcXfJvb1qN5FG6NBoJQXGJ4SxRX9Ynl+WTZ5RWVWl+OV1u8rYsH6PH5zcXfSukRYXY5LaRAo5SUevLI3AX7CAx9t0yeOXazSVsOfPthKXPs2TB/tHReIa9MgUMpLxLZrwx/H9+Lb3YV8tOmA1eV4lZdX5JBVUMojE/sQGuTU41duSYNAKS9y0wXdGJTQnkc+2aFdVbtITmEpzy3NYkK/zozq7dlPEDdEg0ApL+LnJzz+0/MorbDxt093WF2Ox6upMfzx/S0EB/jx8FV9rC6nxWgQKOVlUmLCuf2yJD7adJDlu7RTOme8tiaX9fuO89BVfYiOCLG6nBajQaCUF7p9RE+So9sy44MtFJdVWV2OR8orKuOJL3ZxaUoUPx3kPc8M1EeDQCkvFBzgzzPXD+BoaSV//WS71eV4nJoaw4wPtuAnwj+u7efxnco1RoNAKS/VL74dU0ck8cH3B1i8/bDV5XiU19bksir7GPdP6EVc+zZWl9PiNAiU8mLTRiTRp0sED3y4Ve8iaqLsghL+sWgnI1KjuHGIbwyLq0GglBcLCvDj6ev7c/K0jfs/2KoPmjWi0lbDXW9vIiw4gCeuO8/rTwmd4VQQiEhHEflSRLIcPzvUs06qiGyq9TopInc5lj0sIgdqLZvgTD1KqR/r1TmC+y5PZcmOI7yxdr/V5bi1Z5fuZtuBkzx2bT+iw733LqG6nG0RzACWGmOSgaWO6f9hjMk0xgwwxgwABgNl2AewP+NfZ5YbYxbV3V4p5bxfX9ydS5I78bdPd7D7SInV5biltTnHmPX1Hq5Pj+fyPp2tLqdVORsEE4H5jvfzgasbWX8UsMcYk+vkfpVSzeDnJzx9fX/CQwK4863vdZzjOo6WVnDngu9JjAzjIS9+cKwhzgZBjDHmEIDjZ3Qj698AvFVn3jQR2SIi8+o7tXSGiEwRkQwRySgsLHSuaqV8UHR4CE/9rD+7Dpfw6Gf61PEZNTWGu9/exPGyKv5z4yDaBntfX0KNaTQIROQrEdlWz2tic3YkIkHAT4B3a82eBfQEBgCHgKcb2t4YM9sYk26MSY+KimrOrpVSDiNSo/nt8B68vmY/H36fb3U5bmHWN3tYkXWUh6/q43XdSzdVo9FnjBnd0DIROSIiscaYQyISC5ztefbxwEZjzJFan/3DexF5Gfi0aWUrpc7VfZen2oda/GArvTpH0DvWN7/8AFbvOcrTSzK5qn8XJg3panU5lnH21NBCYLLj/WTg47OsO4k6p4Uc4XHGNcA2J+tRSjUiwN+P528cSERIIL97fQPFp32zC4q8ojKmvrGRHlFt+cc1fX3mVtH6OBsEjwNjRCQLGOOYRkS6iMgPdwCJSKhj+Qd1tn9SRLaKyBZgBHC3k/UopZogOjyEF34xiPzjp7nn7U1U1/jW8wVllTZu/W8GthrD7JsHE+4lg9CfK6euihhjjmG/E6ju/IPAhFrTZUBkPevd7Mz+lVLnLj2xI3/5SR/+/NE2Hlu0kwevTLO6pFZhjOG+d7ew+0gJ8245nx5Rba0uyXK+d3lcKfWDm4d2Y09BKXNW7qVHVFtuvMD7u1T491dZfLb1EDPG9+Ky1MZudPQNGgRK+bgHr+jNvmOneOjjbXSLDOWipE5Wl9RiFqzbz7NLs7hucDy/Hd7D6nLchvY1pJSPC/D34/lJA+kZ1ZbbXtvAtgPFVpfUIpbvKuCBj7YxPCWKx3yga+nm0CBQShEeEsirvzqfiDaB/HLeOrILSq0uyaU25Z3g9jc20js2nBd+MYhAf/3qq03/NJRSAMS2a8Prv7kAP4Ffzl3LgROnrS7JJbbmF/PLuWvpFB7EvFvO98knhxujQaCU+kH3TmHM/9UQSips/OLlNR4fBtsPFnPT3LWEhwTy1q1DfapH0ebQIFBK/Y8+Xdox/1dDOFZayfUvfsf+Y2VWl3ROdh46yU1z1hIa5M9btw4lvkOo1SW5LQ0CpdSPDErowJu3DuVUpY2fvbTa464ZrM05xvUvfUdwgD0EEiI1BM5Gg0ApVa9+8e1YMGUo1TWG61/6jg25RVaX1CSLtx/m5nnriAoP5r3fDSOxU5jVJbk9DQKlVIN6dY7gnd8OIyIkgEkvr+XjTQesLqlBxhjmr97H717fQFpsBO/ddqGeDmoiDQKl1Fn1iGrLh7dfxICu7Zm+YBPPLMl0u76Jyququffdzfxl4XZG9ormzVsvoGNYkNVleQwNAqVUozqEBfH6ry/gusHxPLcsm5vmrOXIyXKrywJg/7Eyrn1hNR9+f4C7R6cw++Z0QoP0FtHm0CBQSjVJUIAfT113Hk/+9Dw25Z1g3L+/ZenOI41v2EKMMby+Jpfxz35L/vEy5k5OZ/roZPz89Inh5tIgUEo1mYhw/fld+eSOi+ncrg2/np/B1Dc3cri4dVsHeUVl3DR3LQ9+tI2BCR1YNP0SRvaKadUavIkY417n+poiPT3dZGRkWF2GUj6tvKqal77J4YWvswnwE6aPTuaXwxIJCfRvsX2eKKtk5vJs5q/OJdBfeOCKNCYN6ar9BjWRiGwwxqT/aL4GgVLKGfuPlfHwJ9tZtquATm2DufWS7vxiaDeXduVwtLSCt9bu5+UVOZRU2LhuUDz3jE0htl0bl+3DF7RIEIjIz4CHgd7AEGNMvd/OIjIOeBbwB+YYY86MZNYReBtIBPYB1xtjjje2Xw0CpdyLMYa1e4uYuTybFVlHiQgJ4Mr+Xbh6QBzp3Tqc03n7quoaMvYdZ8H6/SzaeoiqasPIXtH8YVwqvTr77jjLzmipIOgN1AAvAb+vLwhExB/YjX2oynxgPTDJGLNDRJ4Eiowxj4vIDKCDMeaPje1Xg0Ap97Up7wSvrNrLku1HOF1VTZd2IQztGcngbh0Y2LUDCZGhhAX5/+h0TnFZFVkFJew+UsqqPUf5dnchJeU2woMD+OngeG4a2o2kaB1NzBkNBYGzQ1XudHz42VYbAmQbY3Ic6y4AJgI7HD8vc6w3H/gaaDQIlFLua0DX9jx7w0BOVdj4cscRPt92iG93F/LBxv97GK1NoD9R4cEAVNiqOV1Zzcly2w/Lo8KDGd+3MyN7RXNJchRh2mNoi2qNP904IK/WdD5wgeN9jDHmEIAx5pCINDhunIhMAaYAJCR4/3B6Snm6sOAArh4Yx9UD4zDGsL+ojE15JzhcXE5BSQVHSyvwEyE4wI/gAD/iOrQhKbotSVHhdO3YRi8At6JGg0BEvgI617PoAWPMx03YR31/m80+H2WMmQ3MBvupoeZur5SyjojQLTKMbpHa7487ajQIjDGjndxHPtC11nQ8cNDx/oiIxDpaA7FAgZP7Ukop1Uyt8UDZeiBZRLqLSBBwA7DQsWwhMNnxfjLQlBaGUkopF3IqCETkGhHJB4YBn4nIYsf8LiKyCMAYYwOmAYuBncA7xpjtjo94HBgjIlnY7yp63Jl6lFJKNZ8+UKaUUj6iodtHta8hpZTycRoESinl4zQIlFLKx2kQKKWUj/PIi8UiUgjknuPmnYCjLizHU/jicfviMYNvHrcvHjM0/7i7GWOi6s70yCBwhohk1HfV3Nv54nH74jGDbx63Lx4zuO649dSQUkr5OA0CpZTycb4YBLOtLsAivnjcvnjM4JvH7YvHDC46bp+7RqCUUup/+WKLQCmlVC0aBEop5eN8KghEZJyIZIpItmOMZK8jIl1FZLmI7BSR7SIy3TG/o4h8KSJZjp8drK7V1UTEX0S+F5FPHdO+cMztReQ9Ednl+Dsf5u3HLSJ3O363t4nIWyIS4o3HLCLzRKRARLbVmtfgcYrI/Y7vtkwRubw5+/KZIBARf2AmMB5IAyaJSJq1VbUIG3CvMaY3MBSY6jjOGcBSY0wysNQx7W2mY+/q/AxfOOZngS+MMb2A/tiP32uPW0TigDuBdGNMX8Af+xgn3njMrwLj6syr9zgd/8ZvAPo4tnnB8Z3XJD4TBMAQINsYk2OMqQQWABMtrsnljDGHjDEbHe9LsH8xxGE/1vmO1eYDV1tSYAsRkXjgCmBOrdnefswRwHBgLoAxptIYcwIvP27sIyu2EZEAIBT7iIded8zGmG+BojqzGzrOicACY0yFMWYvkI39O69JfCkI4oC8WtP5jnleS0QSgYHAWiDGGHMI7GEBRFtYWkv4N/AHoKbWPG8/5h5AIfCK45TYHBEJw4uP2xhzAPgnsB84BBQbY5bgxcdcR0PH6dT3my8FgdQzz2vvnRWRtsD7wF3GmJNW19OSRORKoMAYs8HqWlpZADAImGWMGQicwjtOiTTIcU58ItAd6AKEichN1lblFpz6fvOlIMgHutaajsfepPQ6IhKIPQTeMMZ84Jh9RERiHctjgQKr6msBFwE/EZF92E/5jRSR1/HuYwb773S+MWatY/o97MHgzcc9GthrjCk0xlQBHwAX4t3HXFtDx+nU95svBcF6IFlEuotIEPYLKwstrsnlRESwnzPeaYx5ptaihcBkx/vJwMetXVtLMcbcb4yJN8YkYv97XWaMuQkvPmYAY8xhIE9EUh2zRgE78O7j3g8MFZFQx+/6KOzXwbz5mGtr6DgXAjeISLCIdAeSgXVN/lRjjM+8gAnAbmAP8IDV9bTQMV6MvUm4BdjkeE0AIrHfZZDl+NnR6lpb6PgvAz51vPf6YwYGABmOv++PgA7eftzAX4FdwDbgNSDYG48ZeAv7dZAq7P/j//XZjhN4wPHdlgmMb86+tIsJpZTycb50akgppVQ9NAiUUsrHaRAopZSP0yBQSikfp0GglFI+ToNAKaV8nAaBUkr5uP8PYlpddXRrILkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.plot(np.sin(np.linspace(0,2*np.pi,100)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "64f5b287",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from environment import World\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import math\n",
    "import os\n",
    "os.environ[\"KMP_DUPLICATE_LIB_OK\"]=\"TRUE\"\n",
    "from torch.utils.tensorboard import SummaryWriter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8f0c129d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mkdir(path):\n",
    "    folder = os.path.exists(path)\n",
    "\n",
    "    if not folder:\n",
    "        os.makedirs(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "52b0d7b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_location_3d(x_uav, y_uav,z_uav,t, x_user, y_user,z_user, savepath,s):\n",
    "    x_uav = np.transpose(x_uav)\n",
    "    y_uav = np.transpose(y_uav)\n",
    "    h_uav = np.transpose(z_uav)\n",
    "    plt.figure(facecolor='w', figsize=(5, 5))\n",
    "    for index in range(world.urban_world.Build_num):\n",
    "        x1 = world.HeightMapMatrix[index][0]\n",
    "        x2 = world.HeightMapMatrix[index][1]\n",
    "        y1 = world.HeightMapMatrix[index][2]\n",
    "        y2 = world.HeightMapMatrix[index][3]\n",
    "        XList = [x1, x2, x2, x1, x1]\n",
    "        YList = [y1, y1, y2, y2, y1]\n",
    "        plt.plot(XList, YList, color='#708090')\n",
    "    for i in range(user_num):\n",
    "        if s[user_num+i]==1.0:\n",
    "            plt.scatter(x_user[i], y_user[i], c='red', marker='^')\n",
    "        else:\n",
    "            plt.scatter(x_user[i], y_user[i], c='black', marker='^')\n",
    "    for i in range(uav_num):\n",
    "        plt.plot(x_uav[i][0:t+1], y_uav[i][0:t+1], c='blue', marker='.')\n",
    "        plt.plot(x_uav[i][t], y_uav[i][t], c='green', marker='.')\n",
    "    new_ticks = np.linspace(0, 10, 6)\n",
    "    plt.xticks(new_ticks,['0','200','400','600','800','1000'])\n",
    "    plt.yticks(new_ticks,['0','200','400','600','800','1000'])\n",
    "    plt.xlabel('X(m)')\n",
    "    plt.ylabel('Y(m)')\n",
    "    plt.savefig(savepath)\n",
    "    #plt.show()\n",
    "    plt.close()\n",
    "\n",
    "    fig = plt.figure()\n",
    "    ax = fig.gca(projection='3d')\n",
    "    for index in range(world.urban_world.Build_num):\n",
    "        x1 = world.HeightMapMatrix[index][0]\n",
    "        x2 = world.HeightMapMatrix[index][1]\n",
    "        y1 = world.HeightMapMatrix[index][2]\n",
    "        y2 = world.HeightMapMatrix[index][3]\n",
    "        ax.bar3d(x1, y1, 0, world.urban_world.side, world.urban_world.side, world.HeightMapMatrix[index][4]*100, shade=True,\n",
    "                 color='#708090')  #\n",
    "\n",
    "      #  ax.text(x1, y1, world.HeightMapMatrix[index][4]*1000 + 2, str(np.floor(world.HeightMapMatrix[index][4]*1000)))\n",
    "    for i in range(uav_num):\n",
    "        ax.plot3D(x_uav[i][0:t + 1], y_uav[i][0:t + 1], h_uav[i][0:t + 1] * 100,\n",
    "                  c='blue')  # , marker='.', linewidth=3.5, markersize=7.5)\n",
    "        ax.scatter(x_uav[i][t], y_uav[i][t], h_uav[i][t] * 100, c='green', marker='.')  # , markersize=12.5)\n",
    "    for i in range(user_num):\n",
    "        if s[user_num + i] == 1.0:\n",
    "            ax.scatter(x_user[i], y_user[i], z_user[i], c='red', marker='^')\n",
    "        else:\n",
    "            ax.scatter(x_user[i], y_user[i], z_user[i], c='black', marker='^')\n",
    "    new_ticks = np.linspace(0, 10, 6)\n",
    "    plt.xticks(new_ticks,['0','200','400','600','800','1000'])\n",
    "    plt.yticks(new_ticks,['0','200','400','600','800','1000'])\n",
    "    plt.xlabel('X(m)')\n",
    "    plt.ylabel('Y(m)')\n",
    "    ax.set_zlabel('Height(m)')\n",
    "    ax.set_zlim(0, 135)\n",
    "    #plt.show()\n",
    "    plt.savefig(result_path + '3DUrban_' + savepath[8:])\n",
    "    plt.close()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fbe4190a",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Actor(nn.Module):\n",
    "    def __init__(self, state_dim, action_dim, net_width, maxaction,user_num):\n",
    "        super(Actor, self).__init__()\n",
    "        self.fc_t_sprend = nn.Linear(1, user_num)\n",
    "        self.fc_loc_sprend = nn.Linear(2, user_num)\n",
    "        self.l1 = nn.Linear(state_dim - 3 + user_num * 2, net_width)\n",
    "        self.l2 = nn.Linear(net_width, net_width)\n",
    "        self.l3 = nn.Linear(net_width, action_dim)\n",
    "        self.maxaction = maxaction\n",
    "    def forward(self, state):\n",
    "        t = self.fc_t_sprend(state[:, -1:])\n",
    "        loc = self.fc_loc_sprend(state[:, -3:-1])\n",
    "        a = torch.cat([t, loc, state[:, :-3]], 1)\n",
    "        a = torch.tanh(self.l1(a))\n",
    "        a = torch.tanh(self.l2(a))\n",
    "        a = torch.tanh(self.l3(a)) * self.maxaction\n",
    "        return a\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1167f2ef",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'mkdir' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-4b5812ce5f12>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0muav_h\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.95\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mresult_path\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Test_'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0muser_num\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'/'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0mmkdir\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult_path\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[0mT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m200\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[0mLength\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m10\u001b[0m    \u001b[1;31m#1km\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'mkdir' is not defined"
     ]
    }
   ],
   "source": [
    "# test\n",
    "uav_num = 1\n",
    "user_num = 10\n",
    "uav_h = 0.95\n",
    "result_path = 'Test_'+str(user_num)+'/'\n",
    "mkdir(result_path[:-1])\n",
    "T = 200\n",
    "Length = 10    #1km\n",
    "Width = 10\n",
    "R = 5\n",
    "cover_r = 2.5\n",
    "s_dim = uav_num*2 + user_num*2 + 1\n",
    "a_dim = uav_num*2\n",
    "dist_max = 0.5\n",
    "max_action = np.array([math.pi,dist_max/2])\n",
    "\n",
    "test_episode = 25\n",
    "world = World(length=Length, width=Width, uav_num=uav_num, user_num=user_num, t=T, uav_h=uav_h,users_name='Users_'+str(user_num)+'.txt')\n",
    "model = Actor(s_dim, a_dim, 200, max_action,user_num)\n",
    "model.load_state_dict(torch.load('Model_'+str(user_num)+'/td3_actor8000.pth',map_location={'cuda:1':'cuda:0'}))\n",
    "model.eval()\n",
    "\n",
    "coverage_set = np.zeros([T + 1, user_num])\n",
    "energy_set = np.zeros([T + 1, uav_num])\n",
    "\n",
    "success_rate = 0.0\n",
    "Complete_time = np.zeros(test_episode)\n",
    "fly_time = 0\n",
    "hover_time = 0\n",
    "Energy = np.zeros(test_episode)\n",
    "Distance = np.zeros(test_episode)\n",
    "Los = np.zeros(test_episode)\n",
    "Max_cover = np.zeros(test_episode)\n",
    "x0_user = np.zeros( world.user_num)\n",
    "y0_user = np.zeros( world.user_num)\n",
    "z0_user = np.zeros(world.user_num)\n",
    "x0_uav = np.zeros([test_episode, T + 1, world.uav_num])\n",
    "y0_uav = np.zeros([test_episode, T + 1, world.uav_num])\n",
    "z0_uav = np.zeros([test_episode, T + 1, world.uav_num])\n",
    "location = np.zeros([3, world.uav_num])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "949a6e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "def power(v):\n",
    "    P0 = 79.8563\n",
    "    Pi = 88.6279\n",
    "    U_tip = 120\n",
    "    v0 = 4.03\n",
    "    d1 = 0.6\n",
    "    s = 0.05\n",
    "    rho = 1.225\n",
    "    A = 0.503\n",
    "    P_h = P0 * (1+3*v**2/U_tip**2)+Pi*math.sqrt(math.sqrt(1+v**4/(4*v0**4))-v**2/(2*v0**2))+0.5*d1*rho*s*A*v**3\n",
    "    return P_h\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "57b81e54",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'test_episode' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-d9cab7ea3ed6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mfor\u001b[0m \u001b[0mepisode\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtest_episode\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m     \u001b[0ms\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mworld\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreset\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[0mfly_dis\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mmove_time\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mfly_energy\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'test_episode' is not defined"
     ]
    }
   ],
   "source": [
    "for episode in range(test_episode):\n",
    "    s, t = world.reset()\n",
    "    fly_dis = 0.0\n",
    "    move_time = 0.0\n",
    "    fly_energy = 0.0\n",
    "    for i, user in enumerate(world.Users):\n",
    "        x0_user[i] = user.x\n",
    "        y0_user[i] = user.y\n",
    "    for i, uav in enumerate(world.UAVs):\n",
    "        location[0][i] = uav.x\n",
    "        location[1][i] = uav.y\n",
    "        location[2][i] = uav.h\n",
    "    x0_uav[episode][0] = location[0]\n",
    "    y0_uav[episode][0] = location[1]\n",
    "    z0_uav[episode][0] = location[2]\n",
    "    num_fa2 = 0\n",
    "    sum_reward = 0\n",
    "    actions = np.zeros(2 * uav_num)\n",
    "    done = False\n",
    "    while not done:\n",
    "        coverage_set[t] = s[:user_num]\n",
    "        s_shape = torch.unsqueeze(torch.FloatTensor(s), 0)\n",
    "        with torch.no_grad():\n",
    "            actions = model(s_shape)[0].detach().numpy()\n",
    "        fly_dis = fly_dis + (actions[1] + dist_max/2)\n",
    "        fly_energy += power((actions[1] + dist_max/2)*100/2.5)*2.5\n",
    "        s_, r, done, t,terminal = world.step_inside(actions, s, t)\n",
    "        sum_reward += r\n",
    "        move_time += 1.0\n",
    "        s = s_\n",
    "        for i, uav in enumerate(world.UAVs):\n",
    "            location[0][i] = uav.x\n",
    "            location[1][i] = uav.y\n",
    "            location[2][i] = uav.h\n",
    "            #print(uav.x,uav.y,uav.h)\n",
    "        x0_uav[episode][t] = location[0]\n",
    "        y0_uav[episode][t] = location[1]\n",
    "        z0_uav[episode][t] = location[2]\n",
    "\n",
    "        #writer.add_scalar('sum_reward%s'% str(episode),sum_reward, t)\n",
    "        if done:  # and episode % 5 == 0:\n",
    "            print('Episode:', episode, 'Epoch:', t, 'total_Reward: %.4f' % sum_reward, 'reward: ', r, 'fa', world.fa,\n",
    "                  'raw_r', world.r)\n",
    "            print('sum_cover:', world.sum_cover,'hovering_time:',world.hovering_time)\n",
    "            print(world.service_time)\n",
    "            print(world.LoS_time)\n",
    "            print('max_cover_time:',world.max_cover)\n",
    "    #if episode == 6:\n",
    "    draw_location_3d(x0_uav[episode], y0_uav[episode],z0_uav[episode], t, x0_user, y0_user,z0_user,\n",
    "                result_path + 'Test_UAVPath_Users_%s.pdf' % str(episode).zfill(2), s)\n",
    "    if world.terminal:\n",
    "        success_rate += 1.0\n",
    "    Complete_time[episode] = t*2.5 + world.hovering_time\n",
    "    fly_time += t*2.5\n",
    "    hover_time += world.hovering_time\n",
    "    Energy[episode] = fly_energy + world.hovering_time * 168.4842\n",
    "    Distance[episode] = fly_dis\n",
    "    Los[episode] = world.LoS_time\n",
    "    Max_cover[episode] = world.max_cover\n",
    "    print(move_time)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97b986de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i in range(1,T+1):\n",
    "print('#', success_rate)\n",
    "print('##', Complete_time)\n",
    "print('###',Energy)\n",
    "print('#', success_rate / test_episode)\n",
    "print('##', sum(Complete_time)/test_episode)\n",
    "print('###',sum(Energy)/test_episode)\n",
    "print('Distance:', sum(Distance)/test_episode)\n",
    "print('Los time:', sum(Los)/test_episode)\n",
    "print('fly_time: ',fly_time/test_episode,'   hover_time: ',hover_time/test_episode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b825dea",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
